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Shapes From Pixels 

Mitra Fatemi, Arash Amini, Loic Baboulaz, and Martin Vetterli 


Abstract —Continuous-domain visual signals are usually cap¬ 
tured as discrete (digital) images. This operation is not invertible 
in general, in the sense that the continuous-domain signal cannot 
be exactly reconstructed based on the discrete image, unless it 
satisfies certain constraints ( e.g ., bandlimitedness). In this paper, 
we study the problem of recovering shape images with smooth 
boundaries from a set of samples. Thus, the reconstructed image 
is constrained to regenerate the same samples (consistency), as 
well as forming a shape (bilevel) image. We initially formulate the 
reconstruction technique by minimizing the shape perimeter over 
the set of consistent binary shapes. Next, we relax the non-convex 
shape constraint to transform the problem into minimizing the 
total variation over consistent non-negative-valued images. We 
also introduce a requirement (called reducibility) that guarantees 
equivalence between the two problems. We illustrate that the re¬ 
ducibility property effectively sets a requirement on the minimum 
sampling density. One can draw analogy between the reducibility 
property and the so-called restricted isometry property (RIP) in 
compressed sensing which establishes the equivalence of the lo 
minimization with the relaxed A minimization. We also evaluate 
the performance of the relaxed alternative in various numerical 
experiments. 

Index Terms —Binary images, Cheeger sets, measurement- 
consistency, shapes, total variation. 

I. Introduction 

Sampling is at the heart of all digital signal acquisition 
devices. To store and process the data, we need to convert 
continuous-domain signals into a sequence of numbers. Since 
the continuous-domain signal will no longer be available after 
this stage, we expect the sequence to provide an exact or at 
least a fair representation of this signal. The Shannon sampling 
theory and its variations consider sampling strategies for the 
class of signals living in a shift-invariant space. Interested 
readers are referred to 0, EL 0 for a comprehensive study 
of the topic. Here, we focus on the imaging application carried 
out by a digital camera, in which the optical lens and the sensor 
array are responsible for the sampling procedure (Figure [I]). 
The physics of the device imply that the measured sensor 
values—pixel values of the digital image—are samples of a 
filtered visual signal. The impulse response of the involved fil¬ 
ter is called the point spread function (PSF) and is determined 
by the optical system. On the other hand, the sensor array 
controls the number of samples and the sampling resolution. 

In this paper, we study the problem of recovering 
continuous-domain visual signals from their samples. Within 
the broad range of visual signals, we shall focus only on the 
class of shape signals over a fixed background. A shape is 
mathematically described as the characteristic function of a 
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Fig. 1: Sampling system implemented by a digital camera; the 
effect of the optical lens and the sensor array on 2D visual 
signals can be modeled by filtering followed by sampling in 
space. 


union of a few connected subsets. We call these signals shape 
images. We allow for shapes with arbitrary geometries as long 
as they have smooth boundaries. Examples of shape images 
can be found among artworks such as woodcut prints, cutouts 
and lithographs (Figure [2]). 

The shape images have only two different intensity values 
(0 and 1). However, the filtering effect caused by the PSF 
smooths out sharp intensity transitions and results in mea¬ 
surements with varied intensity levels. In short, the sampling 
process projects binary shape images onto gray-scale discrete 
images (Figure [3]). In this paper, we are interested in a 
strategy to recover the original binary image or a binary 
approximation thereof from measurements. In any case, the 
recovered binary image should be able to regenerate the same 
measurements. This requirement assures us that we cannot 
discriminate between the original and the recovered images 
at least at the output of the sampling block. A reconstruction 
of the original image that satisfies this condition is called 
measurement-consistent or consistent for short H, 0, 0. 
The problem of consistent shape image reconstruction appears 
in applications where the aim is to exactly locate or describe 
the objects in a scene; astronomical imaging, quality moni¬ 
toring in manufacturing, biomedical imaging and high-quality 
artwork rendering are a few examples. 

A. Related Works 

Due to the diverse shape geometries and sharp intensity 
transitions on the boundaries, this class of visual signals, 
like many other real-world signals, are neither bandlimitted 
nor belong to a shift invariant subspace. Hence, the classical 
sampling results do not apply here. A similar scenario happens 
for the class of ID signals studied in Q, 0, known as 
signals with finite rate of innovation (FRI). It is shown that the 
discrete samples can lead to perfect signal recovery, although 
the signals are not necessarily bandlimitted. A generalization 
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Fig. 2: Cutout by Henri Matisse (1952) that can be classified 
as a shape image. 


to 2D FRI signals is presented in 0, 0E CD, with the goal 
of recovering convex polygonal shapes from the gray-scale 
pixels. A different approach is devised in flZl by considering 
the boundary curves in a shape image as the zero-level-sets of 
specific 2D FRI signals. Due to the FRI requirements, exact 
recovery relies on the PSF satisfying the so-called Strang-Fix 
condition. Furthermore, the FRI model admits limited shape 
geometries. 

The shape image recovery can also be viewed as fitting 
boundary curves to the interpolated gray-scale image (high- 
resolution version of the measurements). Such methods are 
widely known as segmentation techniques that fit deformable 
curves to gray-scale images, and include active contour al¬ 
gorithms also known as snakes. Based on the curve models, 
they are classified as point snakes CD, geodesic snakes 
cd, ca and parametric snakes ca, ca. In all cases, the 
segmentation algorithm is formulated by minimizing a snake 
energy functional that depends on the gray-scale image and 
the model of the boundary curves. However, it does not take 
the PSF into account ca. As a consequence, the resulting 
binary image is likely to fail the consistency requirements. 


B. Contributions 

In this work, we propose a method to recover a 
measurement-consistent shape image that has continuously 
twice differentiable ( C 2 ) boundary curves. Our approach is 
direct in the sense that it avoids intermediate curve fitting 
steps, and finds the shape with minimum perimeter. We 
formulate the method as an optimization problem constrained 
by the measurements (i.e. pixels), where the functional is the 
continuous-domain total variation (TV). Ideally, we should 
restrict the search domain to binary images. This leads to 
a non-convex problem which is computationally intractable. 
Hence, we consider the convex relaxation in which the search 
is over the set of all non-negative-valued images. Under a 
minimum resolution requirements (see Definition [4] for an 
explicit explanation), we prove that all the solutions to the 
non-convex problem are minimizers of the convex relaxation 
(see Theorem 0. 

The number of constraints in the convex problem equals 
the number of pixels. However, we demonstrate that when 
the resolution requirement is satisfied, the multiple constraints 
can be replaced with a single one formed by a wisely chosen 
linear combination of them. This reduces the problem to an 



(a) (b) 

Fig. 3: Example of a bilevel image and its acquisition: (a) 
a shape image, and (b) pictorial representation of 10 x 10 
measurements, generated with a bilinear B-spline sampling 
kernel. 


equivalent TV minimization problem with a single constraint, 
which is known in the literature as the Generalized Cheeger 
problem CD. a generalized Cheeger set is a shape with mini¬ 
mum perimeter and a fixed weighted integral. This equivalence 
allows us to apply the existing results that show the Cheeger 
solutions are among the minimizers of the relaxed problem 

CD, EE 

The advantage of our method compared to the FRI works 
in (9) -02) is that, we do not constrain the boundary curves 
by any specific model. Instead, we let the sampling kernel and 
the measurement values decide for them. As a result, there is 
less restriction on the achievable shape geometries. Besides, 
the choice of the PSF is arbitrary, and does not need to satisfy 
the Strang-Fix condition. 


C. Organization of the paper 

We explicitly define the problem and the used notations in 
Section [II] We continue by reviewing the concept of Cheeger 
sets and the existing results in Section [111] In Section [TV] 


we present the theoretical results. We employ the primal- 
dual algorithm of KTl for the numerical approximation of 
the solutions to the convex minimization problem in Section 
|Vl This algorithm enables us to study the performance of 
the proposed shape recovery method in Section [IV] through 
numerical experiments. Finally, we conclude the paper in 
Section [VI] 


II. Problem Definition 

We denote by I(x,y) the continuous-domain image with 
domain U = [0, l] 2 . Also, we represent the discrete (or 
measurement) image by D. When D is the output of anmxm- 
pixel digital camera with PSF f(x,y) = f(—x, —y), we can 
relate the m 2 pixels dij , 1 < < m of D to the image 

I(x,y) as 

1 —,xy 

dij = ^2 * d{X,y) I ( x ,y)=(jT,iT) 

= Hx, y )dxd y , (1) 

where T is the sampling period (Figure [TJ. 

In the consistent image recovery problem, we wish to find 
an approximation I of the original image that regenerates the 
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Fig. 4: Perimeter minimization prevents unnecessary details and extra connected components to appear in the shape. In this 
figure, all shapes are consistent with the measurements in (a) but the shapes in (c) and (d) have higher perimeters due to extra 
details on the boundary and additional component, respectively. 


same measurement pixels. This means that the reconstruction 
error I — I between the original image and its approximation 
is in the null space of the imaging process. Equivalently, the 
two images are perceived as identical by the imaging device. 
Let k = (j — l)m + i, 1 < k < m 2 represent the equivalent 
index of dij in the vertical raster scan of D. Also, let 

fk(x,y) = L ~ \k/m\ , |; - ((k mod to) + 1)) 

indicate the sampling kernel in 0 associated with <4. We 
denote by Cq(D; /i,...,/ m 2 ) the set of all non-negative- 
valued images over the domain Cl that are consistent with 

D = [<4] l<k<m 2 > 


C^(D; /i,...,/ m2 ) 

= {/ e BV(Q),I> 0 ; 


( 2 ) 


JJ If k dxdy = d k , l<fc<m 2 j. 


Here, BV(Cl) stands for the set of functions over Cl with 
bounded variation; i.e., all elements of BV(Cl) have well- 
defined and finite total variation values. 

Consistent image recovery is equivalent to finding an ele¬ 
ment of C^(D; /i,..., fm 2 ). In the consistent shape recovery 
problem, we limit the permissible solutions to the shape 
characteristic functions. Let S be a subset of Cl. We denote by 
Xs{x,y) the characteristic function of S over Cl 


Xs(x,y) 


1, if 0 ,y)eS 
0, if (x, y) £ f2 \ <S. 


We call S a shape if it is the union of a finite number of 
connected subsets of Cl. In this case, we call xs a shape image. 

We assume that the sampling kernel in O is always positive 
and has unit C\ norm 



y) dx d y 


= 1 . 


Consequently, the pixels of the discrete image associated with 
a shape characteristic function shall take continuous values in 
the range [0,1]. An example of a shape characteristic function 
and its associated 10 x 10 discrete image is shown in Figure 
[3] We recall that Figure [3b] provides a pictorial representation 
of the 100 pixels in D and it should not be mistaken with a 
piecewise constant approximation of the original image. 

The consistent shape reconstruction problem is equivalent 
to finding a shape image I = xs(x,y) £ C n ( D; /i, ...,/ m 2 ) 


for the set of to 2 pixels 0 < dk < 1 in D. Among all possible 
candidates, we are interested in shape images with minimum 
perimeter. This way we reject shapes with extra connected 
components and excessive boundary details (see Figure [4]). 

Minimum-perimeter consistent shapes are the global mini- 
mizers of the following problem 


inf Per(«S), 

sen, xseBV 

S*t. I XS C fli •••5 fm 2 \ 


OPo) 


where Per (S) is the perimeter of S. Problem is a 

variational non-convex problem and it is prone to having 
many local minima. This makes it very likely that common 
gradient descent methods get trapped in local minima. While 
in problems of this sort, global minimizers are usually all 
reasonable solutions, the local minima can be blatantly false. 
In the next sections, we show that if the discrete image D 
satisfies a resolution requirement defined in Definition [4] the 
minimum-perimeter consistent shapes are the minimizers of a 
convex relaxation of m Furthermore, we conjecture that 
under this condition, there is a unique minimum-perimeter 
consistent shape which is also the unique solution of the 
convex problem. In the experimental section, we present an 
algorithm for the recovery of this solution. 


III. Cheeger Sets 

An image is called consistent with the measurements if it 
complies with all the constraints in Essentially, each pixel 
of the discrete image enters ( fHj] ) as a constraint, resulting in 
an optimization with many constraints. In addition, we are also 
restricting the search domain to bilevel images, which further 
complicates the minimization task. The simplest scenario of 
having only one single pixel (measurement) is a well-studied 
topic known as the Cheeger problem. There is already a 
rich literature regarding the existence, uniqueness properties, 
regularity (smoothness) of the boundary and numerical eval¬ 
uation of such sets for almost arbitrary kernels /. In this 
section, we present a brief review of the Cheeger problem 
and related results upon which we build our general multi¬ 
constraint minimization problem. The details for the latter will 
be discussed in the next section. 

The Cheeger problem can be directly extended to higher 
dimensions; however, for the purpose of image recovery, we 
focus on 2D signals in this paper. Let Cl be a subset of M 2 . The 
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(a) (b) 

Fig. 5: Example of (a) a generalized Cheeger set, for the 
weight kernels g = 1 and / as displayed in (b). 


Cheeger sets of Ll are defined as those 5 C 0 that minimize 
the ratio of the perimeter over the area, 


Per(<S) 

I Is dx d y ' 


(3) 


It is common to represent Per(<S) in terms of the total vari¬ 
ation of the shape image xs- For this purpose we invoke the 
coarea formula that for a positive function u(x, y) : Ll -A M-° 
implies that 


TV(u)= I f \Vu\dxdy = f Per(£/iq /i)) d/i, 

J Jn Jo 


where 


E(u\n) = {(x,y) e 


u(x,y) > y} 


are the level-sets of u(x,y). This immediately indicates that 

Per (S) = TV(xs)- 

We can expand the definition of a Cheeger set by intro¬ 
ducing two weight kernels in the nominator and denominator 
of 0 03. Indeed, a generalized Cheeger set is a shape 
minimizer of 

inf //nfll V X*l d ^ (4) 

scn, X seBV JJ Q fxsdxdy 

Note that for / = g = 1 we obtain the standard Cheeger 
sets. For a simple domain such as a square, the Cheeger set is 
unique and has a certain shape but depending on the choice 
of the weight kernels, generalized Cheeger sets can have very 
diverse shapes. Figure [5] displays an example of a generalized 
Cheeger set for g = 1 and / as in Figure [5b] 

Except for special choices of / and g , the minimizer of <0 
shall not be unique. Furthermore, the minimization in 0 is 
over a non-convex set, which is computationally intractable in 
general. However, an interesting result by Strang in ED (see 
also lED)) shows that all global minimizers of ([4]) (all Cheeger 
sets) are the level-sets of the solution(s) to 

inf //nfl|W|da^ 
ff n fl dx dy 

Note that the search domain in the latter problem consists of 
all non-negative-valued images (not necessarily bilevel), which 
is a convex set. The following statement of this result by (20) 
is more aligned with our approach in the next section. 


Theorem 1 ([20]). Let I be a minimizer of 

inf // q\VI\dxdy , 
ieBVn(f) J 


(5) 


where 

BVn(f) = {/ e BV(Sl), I> 0 


JJ^fldxdy = l). 


Then, for every y > 0 such that the level-set E{I\p) is 
nonempty, 


//E(7 ;/ ,)/d®d2/ XE(7;M) 

is also a minimizer of 0. 


( 6 ) 


In a nutshell, Theorem [T] states that the minimizer set of 0 
is closed under level-set evaluation; i.e., normalized (scaled) 
non-empty level-sets of a minimizer also belong to the set of 
minimizers. This helps in finding a bilevel solution to ([4]), as 
finding any minimizer of the convex problem 0 necessarily 
leads to (at least) a bilevel image. 

Another result proved in (20) indicates that the Cheeger 
sets are closed under set union. This immediately establishes 
the existence of a unique maximal Cheeger set that contains 
all the other ones (23). Thus, we can remove the inherent 
ambiguity caused by non-uniqueness of the solutions to ^ by 
searching for the maximal set. However, finding the maximal 
set is not generally easy by considering the minimization of 
0. A regularization technique is proposed in (23]] that applies 
asymptotically vanishing penalty terms to the cost function 0 
and achieves the maximal set at the limit of the minimizers. 
Based on this idea, a numerical method is introduced in El 
that approximates the maximal Cheeger set on a finite grid. 
The method is robust to discretization as the approximations 
converge point-wise to the continuous-domain Cheeger set 
when the grid resolution increases. 

As a final note, we discuss the influence of the weight 
kernels / and g. In fact, Cheeger sets consist of smooth 
C 2 boundaries, irrespective of the choice of / and g fT9) . 
Nevertheless, it is known that the curvature of the boundaries 
is tightly controlled by these weight kernels. Formally, at each 
boundary point we have that 03 


M < 


J{S) sup / + sup jjVffjj 

inf g 


(7) 


where n stands for the curvature and J{S) is the cost value of 
the Cheeger set determined by the ratio in 0. As we set g = 1 
in the rest of the paper, the effective bound on the curvature 
simplifies to \k\ < J(S) sup/. We will just briefly comment 


on employing a non-constant weight kernel g in Section VI 


IV. Consistent Shape Recovery 


Let us consider the problem <E) for the case where the 
measurement image D = [dk]i<k<m 2 consists of more than 
one pixel. Similar to the single-measurement setting, non¬ 
convexity of the problem is a computational barrier. Therefore, 
we opt to use a convex relaxation in the form of 


inf 

lec n ( D; 


II, 


|V/| cfrch/. 


(Pi) 
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By extending the search domain from binary (bilevel) shapes 
to all non-negative-valued images, the problem becomes con¬ 
vex. Nevertheless, due to multiple measurement constraints, 
this scenario obviously deviates from the conventional Cheeger 
problem. 

In Theorem [2] we show that under certain conditions, the 
minimization in E) constrained by multiple measurements 
can be replaced with a similar minimization subject to a 
single constraint; i.e., we prove that E) could potentially 
have an equivalent Cheeger problem. In fact, we use a wisely 
chosen linear combination of the measurements as the single 
measurement. The interpretation of ( |Pi| ) in form of a Cheeger 
problem automatically implies the existence of a bilevel min- 
imizer ( e.g ., the maximal Cheeger set) for ei Thus, all 
shape minimizers of are also minimizers of the relaxed 
problem (Fj\ . Further, it proves the existence and uniqueness 
of a maximal consistent shape. In Theorem [3j we provide 
a simple test to verify whether a minimizer of E) is the 
maximal shape. This helps us to validate a numerical solution 
obtained via minimizing E) —which might not have a unique 
minimizer—as a binary consistent shape. 

The mathematical requirements for the equivalence of EJ 
with a Cheeger problem (existence of a suitable linear combi¬ 
nation of the measurements) is stated in Definition [4} essen¬ 
tially these requirements imply that the sampling density used 
for obtaining the measurement image needs to be fine enough. 


A. Theoretical Results 


We start by defining the maximal consistent shape. 

Definition 1. A maximal consistent shape with minimum 
perimeter, or MCSMP in short, is a solution to E) whose 
support contains the support of all other minimizers of El 

Note that a MCSMP does not always exist. In general, the 
support union of two minimizers of E) does not necessarily 
generate a minimizer by scaling. It is evident by this fact 
that the claimed equivalent Cheeger problem plays a crucial 
role in our results. In Definition [4] below we will describe the 
sufficient conditions that enable us to associate El» or E) 
to a Cheeger problem. 

Before stating Definition [4| we introduce a few notations 
used in the rest of this section. As we need to linearly combine 
the measurement constraints, we represent the n-dimensional 
coefficient set for the convex combinations by A n : 


A A 


> A n ) e 


0 < A, 


n 

1 >=i}- 




For non-negative-valued images, a zero measurement can 
only happen when the image vanishes over the support of 
the corresponding sampling kernel. Thus, we can exclude the 
support region from our search domain. 

Definition 2. For the measurements D = [dk]i<k<m 2 cor¬ 
responding to the pixels 0 < dk < 1 and sampling kernels 
/i,..., f m 2 , let p denote the number of non-zero pixels and 


A = {i\di > 0} = {ai,... ,a p } 


stand for the index set of active pixels. We define the reduced 
domain by 

= fi r (D;/i,... ,/ m 2 ) = fi\ U supp(fi). 

i e (1,..., m 2 } \ A 

Definition 3. For the measurements D = [dk\i<k<m 2 > sam ~ 
pling kernels /i, ..., f m 2 , and a vector A E A p , we define the 
reduced kernel f x :Q r M-° by 

/ a =(e a ^)/(ew4 

k= 1 k =1 

Here, p, Vt r , and ak are as defined in Definition [2] 

Now, we are prepared to state the Cheeger problem equiv¬ 
alence requirements. 

Definition 4. As before, let D = [dk]i<k<m 2 b e the measure¬ 
ments captured by sampling kernels / 1 ,..., f m 2 with 0 < di, 
leading to p, A, and as in Definition [2] For an arbitrary 
A E A p , we define I x = a\s to be the solution of ([5]) 
corresponding to the maximal Cheeger set S with f x , when the 
domain is restricted to Q r . We call (D; / 1 ,..., f m 2 ) reducible 
if A can be partitioned into K\ and K 2 such that 

(i) Vfc e Ki, X e Ap, Afc =a0 : JJ Qv I x f k dxdy < d k , 

(ii) Vfc € K 2 , A € A p : // Q I x f k dxdy < d k . 

It was explained earlier that the measurements di obtained 
from a binary shape through normalized sampling kernels 
satisfy 0 < di < 1. The requirements in Definition [4] simply 
indicate that the maximal Cheeger solution corresponding to 
any convex combination of the kernels except a given one, 
should result in a strictly smaller measurement observed by the 
excluded kernel. Intuitively, we expect the Cheeger solution to 
have less contribution over the support of the excluded kernel. 
However, there are some exceptions; imagine the case where 
the support of a 3 x 3 block of measurement kernels completely 
coincide with the interior of the binary shape. Thus, we shall 
have a block of all-one measurements. Now, it is likely that the 
maximal Cheeger set corresponding to a linear combination of 
the 8 surrounding kernels (but missing the middle one) using 
symmetric weights fully covers the support of the kernel in 
the middle. Hence, measuring this solution via the middle 
kernel results in di = 1, instead of being strictly less than 
1. The partitions K\ and K 2 in Definition [4] are introduced to 
distinguish between the ordinary ( K \) and exceptional (Kf) 
cases. We postpone further discussion and clarifications about 
this definition to Section HY-B I 

Theorem 2. Let (D; / 1 ,..., f m 2 ) be reducible according to 
Definition [5] Then, all solutions of the non-convex problem 
(Pp) are included in the minimizers of its convex relaxation 
(Pi). Moreover, the solution set of El contains a unique 
MCSMP. 

Our proof of Theorem [2] relies on the following lemma, the 
proof of which is provided in the appendix. 

Lemma 1. For a given dimension n and a set {dk}k =1 C M, 
let K 1 , K 2 be a partition of { 1,..., n}, with the possibility of 
K\ = 0 or K 2 = 0, and let v : A n 1 —M n be a continuous 
function that satisfies 
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(i) V A € A„ : A t -«(A) = ELi^4, 

(ii) Vfc € ifi, A € A n : u fe (A) < <4. 

Afc=0 

(in) V & G iC 2 , A G A n : f/c(A) < 

Then, there exists A* G A n such that v(\*) = [di,..., d n ] T . 

Proof of Theorem [2] The main ingredient of the proof 
is to show that under reducibility condition, E) and E) 
can be associated with a Cheeger problem. To show this, 
first note that Cq(D; /i,...,/ m 2 ) is essentially the same as 
C^ r (D r ; f aiy . • •, /a p ), where p, A, f2 r are defined in Defini¬ 
tion [2] and D r = [dk\keA- In addition, for all A G A p , we 
have that 

Cn r (D r ; / 01 CCn r (l;/ A ) = Wn r (/ A ). 

Therefore, any minimizer of 0 that falls inside 
Cn,.(D r ; f a ) is also a minimizer of Et Besides, if 

E) and 0 have a common minimizer, then, all the solutions 
of (Pi\ shall be among the solutions of {5}. This is indeed, 
what we aim to prove. 

Let I x be the maximal Cheeger set solution of ([5]) on £2 r 
corresponding to the weight kernel f x . Consider the function 



We demonstrate that v(-) satisfies the conditions of Lemma [T] 
The first condition directly follows from 



The reducibility property of (D r ; f ai ,..., f ap ) also establishes 
Conditions ^ and ([m]) of Lemma [T] Consequently, there 
exists A* G A p such that 

Vk( A*) =// I X * fau dxdy = d ak , 1 < k < p. 

J Jn r 

This means that the bilevel maximal Cheeger solution I x , 
which minimizes ([5]), is also consistent with the measurements 
D a- Hence, I x * is also a minimizer of (Pi ] ) as well as (Pp ] ); 
i.e., the three problems (5§ with f x over D r , (Pj\ and (Pp\ 
share a minimizer. This proves the first part of the claim. 

As for the second part, note that all minimizing shapes of 
EJ are Cheeger solutions of 0 . Thus, their support should 
be included in the support of the maximal Cheeger solution 
I x *. In words, I x is a MCSMP. ■ 

Theorem [2] states that under reducibility, the solution set 
of (Pj\ is guaranteed to contain a MCSMP. Although we 
believe that the MCSMP is the unique solution of ( |Pi| ) under 
reducibility, it is yet to be proven. However, we introduce a test 
in Theorem [3] to verify whether an obtained solution to (Pj\ 
is the MCSMP. This test helps us in simulation results, where 
we implement a minimization technique and eventually obtain 
a solution with a numerical precision. First, it is difficult to 
make sure whether the result is precisely bilevel, and second, 
even if it is bilevel, is it the MCSMP? 

Theorem 3. Let (D; /i,..., f m 2 ) with di = 1 for some i be 
reducible (at least one measurement equal to one). If the point 



(C) (d) 


Fig. 6 : Violation of the reducibility criterion due to the 
low sampling density, (a) shows the original binary image 
over a 3 x 3 sampling grid. This generates the measurement 
= 0.0523, when the sampling kernels are the shifts of the 
bilinear B-spline kernel in (b). (c) shows the reduced sampling 
kernel f x corresponding to A = [1,1,1, 0,1,1,1,1,1] / 8 , that 
results in the Cheeger solution (d) with levels 0 and 0.9577 
(reproducing the larger measurement d^ = 0.5247). 


values of a solution to E) never exceed 1 , then, this solution 
is the MCSMP and it is binary (non-zero values are all 1). 

Proof. Let I(x,y) < 1 be a solution to (Pi\ , and let i be 
the index of a measurement equal to 1, i.e., di = 1. By 
comparing I(x,y ) < 1 and di = 1 , we conclude that for 
all (x,y) G supp(/i) we should have I(x,y) = 1 (the kernels 
are normalized). If I is the MCSMP, as it takes the value 1, 
it needs to be binary and the proof is complete. Therefore, let 
us assume the MCSMP to be / 7 ^ /. As previously shown, 
the support of / contains the support of /, which obviously 
contains the support of fi. As / is constant over its support 
and is also consistent with measurement di, we should have 
that I(x,y) = 1 for all (x,y) G supp(/i). Thus, / is binary. 
However, this implies that I never exceeds / at any point, 
while they generate the same set of measurements. In turn, 
this suggests that I cannot be less than / on a set of non-zero 
measure. In other words, I and I are essentially equal at all 
points. ■ 

For recovering a binary shape from discrete measurements, 
we infer the following: when the sampling density is high 
enough to provide the reducibility condition for the measure¬ 
ments, the studied convex problem is potentially able to return 
a consistent binary shape with minimum perimeter. Besides, 
the boundary of the output shall be a C 2 curve. 

B. The Sampling Density Requirement 

Earlier, we claimed that the reducibility condition in Defini¬ 
tion [4] is effectively a requirement on the minimum sampling 
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(c) (d) 


Fig. 7: The measurements of shapes with internal holes never 
satisfy the reducibility requirement, no matter how high is the 
measurement density, unless the original domain is replaced 
with the reduced domain, (a) A binary shape image with an 
internal hole, (b) the corresponding 10 x 10 measurements, 
with ^55 = ^56 = ^65 = (Iqq = 0 , (c) reduced kernel / A , 
with equal contributions from the 20 kernels associated with 
pixels on the borders of the central 6 x 6 sub-grid, and (d) the 
maximal Cheeger solution with levels 0 and 0.9891. 


density. Here, we illustrate this intuition by some examples. 

First, we consider the sampling of the shape in Figure 
| 6 a| over a 3 x 3-pixel grid, employing the bilinear B-spline 


sampling kernel depicted in Figure [ 6 b] This generates the 
measurements 


D 


0.5634 0.0523 0.5750 
0.8996 0.9016 0.8882 
0.5247 0.8817 0.5097 


Particularly, we focus on the d 4 measurement pixel (or d \2 in 
the usual matrix indexing format). It is evident that the value 
of this measurement is considerably lower than its neighboring 
measurement pixels. Intuitively, this sharp transition violates 
the resolution requirement. Now, we check the reducibility 
condition: let us exclude the ^4 pixel and apply equal weights 
for a convex combination of the remaining measurements, 
i.e., A = [1,1,1,0,1,1,1,1, l ]/8 G Ag. Figure 6 c depicts 
the reduced sampling kernel / A , and Figure 6 d shows the 


corresponding maximal Cheeger solution. Although d 4 did not 
contribute in this Cheeger solution, we observe substantial 
leakage over its region from the neighboring pixels. Thus, 
(D, / 1 ,..., /g) is not reducible. Oftentimes, sharp transitions 
between neighboring pixel values indicate lack of sufficient 
density for sampling the boundary curve of the shape (possibly, 
parts with high curvature). Similarly, the reducibility condition 
prevents the value of a pixel dropping substantially below its 
neighbors. 

One of the shortcomings of reformulating as <m> using 


a single reduced kernel f x is that the Cheeger solution never 
admits a hole. Figure [7] provides a pictorial explanation. Here, 
we would like to recover the shape image in Figure [7a] from its 
discrete measurements on a 10 x 10 grid. The hole causes four 
vanishing middle pixels (Figure [7b] ), which make it obvious 
that the shape content is 0 in the middle. We now consider 
a reduced kernel by linearly combining (with equal weights) 
only the 20 kernels associated with the pixels on the perimeter 
of the central 6 x 6 sub-grid (Figure [7c]). As claimed, the 
Cheeger solution to ([5} depicted in Figure |7d] has no holes and 
completely covers the middle part. This seems to violate the 
reducibility condition, no matter how high we set the sampling 
density. However, note that we remove the 0 pixels from the 
domain in Definition [4] Therefore, the Cheeger solution over 
the reduced domain is forced to have a hole, although it is not 
considered as hole with respect to the reduced domain. 


The reducibility condition in Definition [4] is a useful guar¬ 
antee for recovering a shape image. However, verifying it 
for a given set of measurements and sampling kernels is a 
combinatorial problem in general. For the purpose of illustra¬ 
tion, we investigate the simple case with 2 x 2 measurement 

pixels. Let D = with elements in [0,1] represent the 

[(22 U4J 

measurement matrix. Without loss of generality, we assume 
that <^4 = p < 1 is the largest element. To verify the 
reducibility condition, we need to exclude each pixel, apply 
an arbitrary convex combination on the rest and check an 
inequality. As we can categorize d 4 to the K 2 set in Definition 
[4] the inequalities when d/± is excluded are trivial. To verify 
other inequalities, note that we can scale all measurements 
by the factor ^ (or any other positive real). In fact, the 
scaling does not affect the support set of the Cheeger solutions. 
Consequently, the reducibility condition for D, boils down to 
a set of inequalities on each °f —, y, y in terms of the other 
two: 


, ^ r,d 2 d 3 
d\ > pZ{ —, —), 

P P 

d 2 >P Y( 

P P 

d 3 > pY( — , —). 
P P 


The symmetries of the problem indicate that the lower-bounds 
on ^2 and d% can be represented using the same function 
(F(-, •)), and the lower-bound Z(-, •) on d\ is symmetric with 
respect to the two inputs. In Figures [ 8 ] and [9] we depict the 
functions Y, Z for two choices of the sampling kernel, namely, 
the box-spline (Figure [ 8 ]) with non-overlapping kernels and 
bilinear B-spline kernels with 50% overlap (Figure [9]). The 
overlap introduces correlation among the neighboring pixels, 
which naturally leads to tighter regions for validity of the 
reducibility condition. This is indicated by larger Y and Z val¬ 


ues. For instance, the measurement set D = 


0.576 

0.216 


0.72 

0.216 


is reducible under the box-spline sampling kernels, but not 
under the bilinear B-spline kernels. 
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(a) (b) (c) 

Fig. 8: The approximate functions (b) Y and (c) Z for the characterization of reducible measurements Y> 2 x 2 when the kernels 
/i, / 2 , fs, are shifts of the box-spline kernel in (a) centered at point (0.25,0.75), (0.25,0.25), (0.75,0.75) and (0.75,0.25), 
respectively. 



(a) (b) (c) 

Fig. 9: The approximate functions (b) Y and (c) Z for the characterization of reducible measurements Y> 2 x 2 when the 
kernels /i, / 2 , / 3 , /a are shifts of the bilinear B-spline kernel in (a) centered at point (0.25,0.75), (0.25, 0.25), (0.75, 0.75) and 
(0.75,0.25), respectively. 


V. Numerical Experiments where 


In this section, we aim at numerically calculating the opti¬ 
mal solution(s) of the convex problem ([Pi}. For this purpose, 
we restrict the simulations to the discrete setting. Below, we 
first explain the equivalent problem in the discrete domain and 
then, present the simulation results. 

A. Discrete Formulation 

For conducting computer simulations, we are limited to 
discrete scenarios. Therefore, we discretize the domain Q (and 
subsequently all the functions defined on Q) with a finite step- 
size h ~ -k for some large integer TV. This will approximate 
Fl and the continuous-domain objects /, /i ,...,/ m 2 by their 
pseudo samples at the 2D grid 


(Y7J) 1 = f li+hj Ii,j lf^<TV, (Q\ 

V Vi A,j | o if i = TV, w 

(V/) 2 ■ = l IiJ+1 - IiJ lf J < N ’ (10) 

y yd h,j if j = n. K ’ 

It is shown that in the asymptotic regime of TV -A oo, the 
results obtained with the discretized model converge to their 
continuous-domain counterpart introduced in ( fPj} l25l . 

One of the standard approaches for solving the associated 
discrete minimization is the gradient descent algorithm, which 
is rather slow in high dimensions (small step-size h). In this 
paper, we use the recently proposed primal-dual algorithm in 
(ZD, EH, which is significantly faster and enjoys convergence 
guarantees. This scheme is based on the weak formulation of 
the total variation: 


{(ih,jh); i,j = 1,2, 

resulting in K JVxiv matrices. In the discretized version, we 
approximate the gradient operator by evaluating the forward 
differences; for instance we approximate V/ with an R NxNx2 
tensor defined as 


(yi)i,j,k = (yi)ij 


( 8 ) 


min 

iec 



= min 
iec 


max 

|C|2<3 


{//-ZdK} 


= min max (— Ldb/C) 
iec |c| 2 <s 


(ID 


with the dual variable ( : D —>> M 2 . Here, div£ stands for 
the divergence and is defined as the negative of the gradient 
adjoint. Each iteration of the optimization algorithm alternates 
between a gradient descent and a gradient ascent on the primal 














9 



(a) (b) (c) (d) 

Fig. 10: The performance of algorithm E) in shape recovery: (a), (b), (c) and (d) show the original shape image displayed 
with the resolution 2000 x 2000 and its approximations using 40 x 40, 50 x 50 and 80 x 80 measurements, respectively. Note 
that the reconstructed images are binary only when the number of measurements is large enough. 


and dual variables, respectively. In short, the update equations 
are as follows: 

c (fc+i) _ Proj B(g) (C fe + <TfcV7 (fc) ), (12) 

I (k+ 1) = Proj Cs2(D;/l) ..., /m2 )(/ (ft) + r fe divC (ft+1) ), (13) 

9k = TF=F=. i T k +1 = OkTk, (Jk+l = 0'k/9k, (14) 

V1 + 4r fe 

i {k+1) = i (fc+1) + e k (i {k+1) - / (fc) ), (15) 


where k represents the iteration index. Here, the notation 
Proj A (-) stands for the orthogonal projection of the argument 
onto the set A and B(g) represents the ball with radius g in 
the space of N x N x 2 tensors: 


B(g) = |u e R NxNx2 


u h, i 


U i,j, 2 — 9i,j 


Hence, Projg^(-) scales only the points outside the 
ball 13(g). The somewhat more complicated projection of 
Proj Cn ( D; j i5 j 2 ^(-) is also implemented using the POCS 
algorithm EZI The initial values I and are arbitrary, 
with 7° = I and time steps Tocro||V || 2 < 1 ID]]. By analogy 
(continuous setting), the divergence in (13] ) shall be the negated 
adjoint of the discrete gradient used in ( fT 2 ] ). For the forward 
difference gradient in equations ([ 8 ])-(^Q|), this leads to 


(divC)zj 


c 1 - r 1 

s i,j 

if 1 < i < TV, 

C7 

if i = 1, 

-clu 

^5 

<S>. 

II 

5-0. 

1 

to 

5-0. 

1 

if 1 < j < N, 

c 2 - 

if j = 1, 

-r 2 

1 

if j = N. 


B. Simulation Results 

In the first experiment, we study the effect of the number 
of measurements on the reconstructed images obtained with 
the proposed algorithm. Recalling the result of the previous 
section, we expect the solution of E) to be binary, given 
adequate number of measurement pixels. In this experiment, 
we employ a shape image with a parametric description, 


composed of a semicircle laid on one side of an equilateral 
triangle (Figure |10a| ). This enables us to precisely access and 
display the image at arbitrary fine resolutions as a reference. 
Figure |10a| shows the image at the resolution 2000 x 2000. 
Figs. |10b| |10c| and |10d| show the solutions of algorithm 


E) with the same resolution applied to the measurement 
of sizes 40 x 40, 50 x 50 and 80 x 80, respectively. All 
measurements are generated with a box-spline kernel. The 
original shape has non-smooth details around the comers and 
thus, to facilitate comparison, we enlarged the reconstructed 
images around these areas. The results reveal that with lack 
of enough measurements, the reconstructed images have more 
than two levels. It seems that the 80 x 80-pixel image provides 
enough measurements to have a binary optimal solution for 

El 



Fig. 11: The average MSE of the recovered binary images 
constrained bymxm samples of random circles with fixed 
radius for two radii r = 0.3 and r = 0.4. 

In a similar experiment, we examine the performance of our 
algorithm in recovering circles from different number of mea¬ 
surements. For this purpose, we run a Monte Carlo experiment 
by generating 20 circles with fixed radius and random centers 
in the image plane. We then consider outputs of the algorithm 
at resolution 600 x 600 constrained with m x m analytic 
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(e) 


(f) 


(g) 


(h) 


Fig. 12: Consistent shape reconstruction with the proposed algorithm: Figs. 12a and 12c display two shapes at the resolution 
1000 x 1000 that will be approximated from 200 x 200 discrete images, generated with biquadratic B-spline sampling kernels. 
Figs. |12b| |12e| and |12f| show the same enlarged sections of the original shape |12a| the corresponding discrete image and 

and 


reconstructed shape, respectively. Figs. 


12d 


12 g 


12h 


display the same for the shape in Figure 


12c 


TABLE I: Quantitative evaluation of the proposed algorithm 
(numbers in dB) 



shape in Figure |l2a| 

shape in Figure |l 2c | 


image 

PSNR 

measurement 

PSNR 

image 

PSNR 

measurement 

PSNR 

proposed 

solution 

29.1507 

58.3316 

43.8839 

63.3429 

linear 

interpolation 

26.4397 

41.4224 

33.7172 

49.1746 


measurements of the circles with box-spline PSFs and different 
values of m. Figure [IT] shows the average mean squared errors 
of the reconstructed images (after thresholding at level 0.5) 
versus m for two different radii. The plots in this figure 
clearly indicate that the algorithm always perfectly recovers 
the circles from mxm measurements when m is greater than 
10. Next, we examine the solutions of <\P\ | ) to 200 x 200- 


pixel discrete images of the shapes depicted in Figures 12a 


and |12c| at resolution 1000 x 1000 (Figure |12c| is taken from 
the middle part of Figure [2]). The sampling kernel for this 
experiment is the biquadratic B-spline. Figure [represents the 
same enlarged sections of the original shapes, the discrete 
images (just for a visual comparison) and their reconstructions 
with the proposed algorithm. The figures demonstrate that both 
reconstructed images are almost binary. Also, Table [T] shows 
the quantitative evaluation of the reconstructed images. In this 
table, we also compare our results with the ones obtained by 
the interpolation of the measurement images with the bilinear 
B-spline kernel, followed by a thresholding at level 0.5. For a 
fair comparison, we also threshold our results to calculate the 
PSNRs. The numbers in this table clearly indicate the success 
of our proposed algorithm for consistent shape reconstruction. 





(b) 


(c) 


Fig. 13: Consistent shape reconstruction with the proposed 


algorithm from highly blurred discrete images: Figure |13a 
shows a 200 x 200 discrete image corresponding to the shape 
in Figure |12a| when the sampling kernels are shifts of an 
stretched biquadratic B-spline with an effective support of 


40 x 40 pixels. Figs. |13b| and |13c| show the same enlarged 
sections (as in Figure |12b| ) of the measurement image and re¬ 
constructed shape, respectively. The recovered image (without 
any thresholding) has a PSNR of 33.8096 dB with respect to 
the original shape and a measurement PSNR of 75.0489 dB. 
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For a given shape image, the resolution requirement in Def¬ 
inition [4] mainly depends on the sampling grid, rather than the 
PSF. To examine this fact, we repeat the experiment in Figure 
12] by regenerating a 200 x 200 discrete image from Figure 
12 a| using a stretched biquadratic B-spline sampling kernel 
with an effective support of 40 x 40 pixels. The result is the 
highly blurred image in Figure |13a| Also, Figure |13b| shows 
the enlarged section equivalent to Figure |12b| The quality 
of the reconstructed image in Figure |13c| (PSNR = 33.8dB) 
confirms that the sampling grid outweighs the choice of the 
PSF in determining the performance. 

Finally in the last experiment, we study the performance 
of the proposed method in a setting severely deficient in the 
number of measurements. For this purpose, we consider a 
recent image by the New Horizons spacecraft in July 2015 
from a moon of Pluto named Hydra. Figure |14a| depicts the 
received measurements. Although a high resolution imager 
is used, due to the long distance of the spacecraft to Hydra 
compared to the size of Hydra, we observe a highly pixelated 
image. According to the available data, the effective PSF 
width of the imager is around 1.5 pixels, which we model by 
a dilated biquadratic B-spline. Figure 14b shows the output 
of the convex program to the measurements by applying 
the approximate PSF. As the measurements are too few, the 
reconstructed image is not bilevel (indeed, it is not unlikely to 
assume the image of Hydra being binary from this distance). 
Nevertheless, it is interesting to note that the obtained multi¬ 
level image is not far from the processed image released by 
NASA in Figure |14c| 


VI. Conclusion 

In this paper, we studied the problem of reconstructing 
a continuous-domain shape image from the samples in a 
gray-scale discrete image. This is essentially equivalent to 
the interpolation of pixels in a way that generates a binary 
image. We formulated this problem as a minimization problem 
where the functional is the continuous-domain total variation 
and the constraints encode the sampling relation between the 
continuous-domain image and the pixels of the discrete image. 
When the search is over binary images, the minimizers will 
be shapes with minimum perimeter and smooth boundaries 
that satisfy the measurements. However, the search over shape 
images is computationally intractable. We introduced the re- 
ducibility condition on the samples of the discrete image and 
proved that when it is satisfied, extending the search domain 
to the non-negative-valued images would not omit any of the 
binary minimizers. The reducibility condition essentially calls 
for smooth changes in the values of the neighboring pixels. 
From this perspective, this is an intuitive requirement on the 
minimum sampling density that is needed for tracking local 
changes in the shape boundaries. 

We conjecture that under the reducibility condition, the 
convex problem has a unique binary solution. Nevertheless, 
we introduced a test to verify whether an obtained solution to 
the convex minimization problem is binary. This test is mainly 
useful in the numerical calculation of the minimizers where 
the recovered solutions might not be precisely bilevel due to 
the numerical precision. 


Our approach in this paper was based on minimization of 
the total variation, but all the results remain valid if we use a 
weighted total variation. A carefully designed weighting kernel 
might locally adjust the recovered shapes and lead to shapes 
with higher mean curvature. 

Appendix 
Proof of Lemma[1] 

In the following, we reserve the notation e f for the canonical 
basis of M n : 


e” = [o,...,o, 1 ,o,...,o] T el". 

ith 

We prove the lemma by induction on n. We set the basis of 
the induction on n = 1. It is trivial to check that Condition ^ 
for n = 1 implies the claim in this case. Next, by assuming 
the validity of Lemma |T] for some n > 1 , we demonstrate the 
validity for n + 1 . 

For the case K[ n ^ = 0, it is not difficult to see that A = 
..., ^q-J T satisfies the requirement. Here, Condition (jij) 
implies that all the inequalities of Condition ([m]) are in fact 
equalities. Hence, we focus on K[ n ^ 7 ^ 0. Without loss of 
generality, we assume that n + 1 G iC{; n+1 \ Next, we will 
try to reduce the (n + 1 )-dimensional problem into a similar 
n-dimensional one with K[ n ^ = K[ n+1 ^ \{n + l} and = 

rAn+l) 

According to Condition (ml), at A = e£ jj; e A n+1 we have 
that 

v;ei^ n+1) \{n + 1 }: Vi (e n n X\) < d % . 

If Ff^ n+1 ) \ {n +1} = 0, set e = Otherwise, set 0 < e < \ 
such that for all i G K[ n ^ \ {n + 1} and all A G A n+ i with 
||A — e^j|| < e (i.e., e-neighborhood of inside A n+ i), 
we have that vi{ A) < d{. The existence of such e follows 
from the continuity of v (and consequently t^s). Furthermore, 
Condition |ih]) implies 14 (A) < di for all i G K ^ n+1 ) and the 
same set of A vectors. In summary, we conclude the existence 
of 0 < e < \ such that 

VI < i < n, A G A n+ i, ||A - e™j{|| < e : u;(A) < d { . 

By taking Condition ^ into account, we observe that 

V A G A n+ i, || A - eJJliH < e : v n+1 (\) > d n+1 . (16) 

In words, the value of u n+ i in a neighborhood of never 
drops below the desired value d n +\. In contrast, the values 
of u n+ i on the facet of the simplex A n+ i opposite to 
(A G A n+ i, A n+ i = 0) are strictly below d n +\ according 
to Condition ([n]). Since u n+ i is continuous, by starting from 
any point on this facet and gradually moving towards 
on the line connecting the two points, v n+ i will eventually 
attain the value d n + 1 . By considering the points on all such 
lines that v n+ i attains the value d n+ 1 for the first time (when 
moving away from the facet towards the vertex {), we shall 
have a manifold intersecting with all the facets except possibly 
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(a) (b) (c) 

Fig. 14: Performance of the proposed method in a setting with limited measurements: the measurement image in (a) is taken 
by the imager of the New Horizons space probe from Hydra. Due to the deficiency of measurements, our reconstruction in 
(b) is not bilevel; yet it is a good match to the processed image (c) released by NASA. 


the studied one. To mathematically represent this manifold we Also, let t G A n be such that ti = 0 for some 1 < i < n. 
employ the following definition: Recalling the definition of u , we have that 


V t G A n : 


P{ t) - inf |/3 G [0,1] 


V7, /3 < 7 < 1 : 


V n +i(lh, ■ • •, 7 tn, 1 - 7 ) < ++i}- (17) 

It is not difficult to apply the continuity of v n +\ to conclude 
the continuity of /3( t) and the fact that 

Vt 6 A n : v n+1 (/3(t)ti, ..., /3(t )t n , 1 -£(t)) — dn+1' 

(18) 


Moreover, we invoke © to demonstrate that /3(t) > i. e ., 

/3( t) is strictly positive for all t G A n . 

Now we are ready to reduce the dimension to n. For this 
purpose, we define the function u : A n h+ M n as 


II 

> 

,. . . , t n ] G A n 



(19) 


«i(/?(t)ii, ... 

, /3(t )t n , 1 - /3(t)) 


’ «i(t) 

u{ t) = 

v n (f(t)ti, ••• 

, /3(t )t n , 1 - /3(t)) 




«t( t) = Vi(\), 


where 


A = 


P(t)h, P{t)t n , 1 - Pit) 


n+1 


K —/ 3 (t) U +1 — /?(t) — 1 => A G A: 


n+l • 


7=1 


7=1 


As ti = 0 results in = 0, the Conditions ^ and ffi]) directly 
carry over to the functions U{ with K[ n) = K[ n+tJ \ {n + l} 
and K^ n) = K^ +1) . 

To sum up, u is a continuous function that satisfies Con¬ 
ditions (p])-([iii|). Therefore, we conclude by the assumption of 
the induction that there exists t* G A n such that 


u{ t*) = [di,...,d n ] T . 


Finally, by plugging this result into (\9\ and using ( p~8] ), we 
. obtain that 


*)+ ...,/3(t*)+ l-/3(t*)) = [d 1 ,...,d n+1 ] T . 


The continuity of ix(t) directly follows from the continuity of 
v and (3. To verify Condition ^ for u note that 
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